Detection of companion galaxies around hot dust-obscured hyper-luminous galaxy W0410-0913

The phase transition between galaxies and quasars is often identified with the rare population of hyper-luminous, hot dust-obscured galaxies. Galaxy formation models predict these systems to grow via mergers, that can deliver large amounts of gas toward their centers, induce intense bursts of star formation and feed their supermassive black holes. Here we report the detection of 24 galaxies emitting Lyman-α emission on projected physical scales of about 400 kpc around the hyper-luminous hot dust-obscured galaxy W0410-0913, at redshift z = 3.631, using Very Large Telescope observations. While this indicates that W0410-0913 evolves in a very dense environment, we do not find clear signs of mergers that could sustain its growth. Data suggest that if mergers occurred, as models expect, these would involve less massive satellites, with only a moderate impact on the internal interstellar medium of W0410-0913, which is sustained by a rotationally-supported fast-rotating molecular disk, as Atacama Large Millimeter Array observations suggest.

called blowout phase 6,7 ), leading to an evolutionary transition from dusty starburst galaxies to ultraviolet (UV)-bright unobscured quasars, and finally to giant elliptical galaxies 6 .
Models predict that the growth of these massive star-bursting objects is driven by galaxy mergers 8,9 , which can deliver a significant amount of gas toward their centers, boosting their star formation rate (SFR) up to thousands of M ⊙ yr −110 , and feeding the growth of their supermassive black holes (SMBH), thereby powering the AGN activity 8 . However, while cumulative number counts in wide-field surveys suggest that Hot DOGs tend to lie in dense large-scale environments 11,12 , direct observations of the mechanisms responsible for their rapid growth are still sparse or only restricted to scales of a few tens of kpc 13 . On the other hand, recent works [14][15][16] have demonstrated the potential of spectroscopic imaging studies across the environment of luminous quasars and protocluster galaxies, finding companion galaxies out to hundreds of kpc.
In this work, we report rest-frame UV and rest-frame far-infrared (FIR) spectroscopic imaging observations of the Hot DOG WISE J041010.60-0913056.2 (hereafter W0410-0913), obtained with the integral field spectrograph Multi Unit Spectroscopic Explorer (MUSE) at the Very Large Telescope (VLT), and with the Atacama Large Millimeter/Sub-millimeter Array (ALMA), respectively. W0410-0913 was originally selected from a WISE all-sky survey 1 and, with a total FIR luminosity of L FIR~2 × 10 14 L ⊙ 17 , and with large stellar 18 (M star ) and molecular gas 19 (M gas ) masses of ≳10 11 M ⊙ , is one of the brightest, most massive and gas-rich galaxies in the distant Universe 19,20 . Extensive restframe FIR observations, tracing several molecular lines, have shown W0410-0913 to be located at a systemic redshift of z = 3.631 18,21 (see details on the redshift determination in the subsection ALMA data in Methods), about only 1.7 Gyr after the Big Bang. Through the combined analysis of the Lyman-α transition of atomic hydrogen (at a restframe wavelength of 1215.67 Å) and the rotational transitions of carbon monoxide, these observations allow us to study both the circumgalactic environment surrounding the hyper-luminous galaxy, on scales of ≳100 kpc, as well as the morphological and kinematic structure of its ISM, on scales of a few kpc, providing information on the multi-scale processes that drive its evolution.

Results
MUSE observations of the circumgalactic environment of W0410-0913. We have obtained deep MUSE observations of W0410-0913 (exploiting a total telescope time of 2 h) as part of an observational campaign aimed at obtaining spatially resolved rest-frame UV spectroscopy observations of a Hot DOG and its halo gas, also called circumgalactic medium (CGM). Due to its wide field of view (FOV) and high sensitivity, MUSE is a suitable instrument to explore the surrounding environment of luminous sources in the distant Universe, and it has carried out observations of both morphology and kinematics of the gas distribution around luminous distant quasars 22,23 through direct imaging of the Lyman-α emission, which is mainly produced through recombination radiation following photo-ionization powered by bright sources 24 .
The spectrum of W0410-0913 observed with MUSE shows highly attenuated rest-frame UV emission. We detect very broad (full widths at half maxima of about FWHM~2800 km s −1 ) blueshifted components of NVλ 1240 and C IVλλ 1548,1550, both offset by more than 2000 km s −1 with respect to the systemic redshift, likely tracing fast AGN-driven outflows 5,13,25 . Also, we find a narrow (FWHM~400 km s −1 ) Lyman-α line, which, similarly to the narrow component of the optical [OIII] λ 5007 line 5 , is blueshifted by about 1400 km s −1 . We find evidence for a small Lyman-α nebula associated to this line, with a maximum projected linear size of about 30 kpc. This is at odds with recent observations of several Lyman-α nebulae showing sizes up to hundreds of kpc detected with MUSE around UV-bright quasars, with typical L bol even lower than W0410-0913 and shorter exposure times 22,23 . Our finding is rather to be ascribed to the large amount of dust in W0410-0913, which prevents Lyman-α and UV photons to escape the galaxy and illuminate the CGM, suggesting that extreme physical conditions of the ISM can have a profound impact on the properties of the CGM accessible through Lyman-α.
The lack of extended diffuse Lyman-α emission dominating the field at wavelengths close to the systemic redshift of W0410-0913 makes it possible to efficiently search for Lyman-α-emitting companion galaxies in the MUSE datacube, and thus provides well-suited conditions to explore the surrounding environment of a Hot DOG on circumgalactic scales, informing on the physics that drives the rapid evolution of this rare population of objects.
Our MUSE observations (reaching a 2σ surface brightness (SB) limit of 5 × 10 À19 erg s À1 cm À2 arcsec À2 in a 1 arcsec 2 aperture and in a single spectral channel) reveal a high overdensity of 24 Lyman-αemitting satellite galaxies associated with W0410-0913 within the MUSE FOV, in a region of about 400 projected physical kpc (diameterscale) centered around the central massive luminous object (see Fig. 1). A compilation of their Lyman-α flux maps is shown in Fig. 1a, providing information on their projected spatial distribution in the field. The adopted techniques and the analysis performed to search for these galaxies (as well as a statistical assessment of the reliability of their detections) are discussed in the subsections describing the MUSE data in Methods.
Archival images from the Hubble Space Telescope (HST) obtained with the WFC3 F160W filter (probing rest-frame emission at around 0.3 μm) show that, despite the modest exposure (only 2400 s), ten among these sources (i.e., approximately 40% of the total) have rest-frame optical counterparts detected with a signal-to-noise ratio (SNR) higher than 4σ, and corresponding SFRs ranging between ≳12−100M ⊙ yr −1 (computed using standard L UV/optical −SFR calibrations; see subsection Measure of the SFR from HST in Methods). Most of the remaining objects show only tentative individual detections at SNR of 2−3σ, but from a stacking analysis, we can infer an average SFR of about 3 M ⊙ yr −1 .
MUSE data show that the companion galaxies discovered around W0410-0913 have offset line-of-sight velocities, traced by the peak of Lyman-α in their spectra (see Methods, where resonant scattering effects are also discussed), in the broad range of about [−6.000; +10.000] km s −1 with respect to the systemic redshift of the central galaxy (see Fig. 1b, c), traced by its molecular gas. 75% of them (18 out of 24) show offsets in a smaller velocity range of about [−2.500; +1.500] km s −1 (see subsection The circumgalactic environment of W0410-0913 observed with MUSE in Methods; see Fig. 1c), suggesting that these satellites are experiencing a strong attraction by W0410-0913.
To quantify the physical significance of the revealed overdensity we compare the Lyman-α luminosity-weighted number density of the companion galaxies associated with W0410-0913, n LyE , with the expected value for field galaxies according to the Lyman-α luminosity functions measured in wide and deep surveys with MUSE 26,27 (see subsection The physical significance of the overdensity in Methods for details on the computation and on the comparison sample). We find that the observed n LyE is a factor of ρ CGM~1 4 higher than the expected value (〈n LyE 〉; see Fig. 2). This is higher than other reported observations of overdensities previously discovered through Lyman-α in the redshift range z2 −7 28-31 , including proto-clusters around radio-galaxies 32,33 and multiple-quasar objects 34 , although a direct comparison with these systems is challenged by the differences in selection criteria, cosmological volumes, and sensitivity limits. Also, we note that the measured overdensity is consistent with cosmological predictions in the proximity of massive, luminous objects like W0410-0913, at the redshift of interest (see subsection The physical significance of the overdensity in Methods).
Altogether, our findings suggest that the circumgalactic environment surrounding W0410-0913 (probed by MUSE on projected physical scales of approximately 400 kpc) is one of the densest regions of the observed Universe. The number of gravitational encounters is expected to be enhanced in such a dense environment 6,17,35,36 , suggesting that galaxy mergers could be an efficient channel for the mass build-up of the massive luminous galaxy.
However, our MUSE observations cannot firmly confirm this scenario, since we do not detect direct unequivocal signs of mergers. A tentative indication of ongoing dynamical interactions is provided by the morphological and kinematic distribution of gas along the line-ofsight as traced by the Lyman-α emission in and around the companion galaxies, that we study by using commonly adopted flux-extraction techniques that enable one to recover flux of extended sources allowing for local kinematic variations 37 . Specifically, we find that the three-dimensional (3D) distribution of circumgalactic gas illuminated by the Lyman-α-emitters dynamically associated with the central luminous galaxy (i.e., with an offset velocity range of about [−2.500; +1.500] km s −1 ) is preferentially distributed over bridge-like shapes on scales up to ≳100 kpc, some of which appear to be aligned along the common axis with the central massive object (see Fig. 3). While these filamentary structures might be induced by the strong gravitational influence of W0410-0913, projection effects combined with the resonant nature of Lyman-α challenge this interpretation.
Further studies of the morphology and kinematics of the distribution of gas in the halo, tracing e.g., bright rest-frame FIR lines (like [CII]) with ALMA and rest-frame optical/near-IR lines with JWST, are needed to discern whether W0410-0913 is in the process of accreting its neighbors, similarly to what is done in highly interacting dense groups of massive galaxies in the local Universe through maps of atomic hydrogen filaments 38 .
Furthermore, the merging scenario is at odds with high-resolution spectroscopic imaging data of the internal gas structure in the ISM of the hyper-luminous central object. Lyman-α-emitting companion galaxies revealed by MUSE in the field of W0410-0913. We applied a spatial Gaussian filtering of 1 pixel (0.2″) to the cutouts and we set a lower bound of SB = 5 × 10 À18 erg s À1 cm À2 arcsec À2 in the color-scale to improve the visualization. The Lyman-α emitters are distributed on scales of about 400 kpc around the W0410-0913 (see the black line, which corresponds to 300 kpc). The size of each cutout varies between 2″ and 4″, and is chosen according to the extensions of the emitting sources. See Table 1 for the definition of IDs and a summary of the measured properties of the Lyman-α emitters. b We color the cutouts shown in panel a according to the offset velocities of the companion galaxies with respect to the systemic redshift of W0410-0913. The redshift of each emitter is traced by the center wavelength obtained through a Gaussian fit of the Lyman-α spectrum. For W0410-0913, at the center, we show the velocity as traced by its Lyman-α, which is blueshifted with respect to the systemic redshift. The offset velocities range between −6.000 km s −1 and 10.000 km s −1 . A histogram of their distribution is shown in panel c. Most of the companion galaxies (18 out of 24) show offsets within a smaller velocity range of about [−2.500; +1.500] km s −1 . These are more likely subject to the gravitational effect of the central massive object. We quantify the significance of the overdensity of companion galaxies discovered around W0410-0913, both in the whole volume defined by the offset velocity range [−6.000; +10.000] km s −1 (green triangle) and in the denser volume defined by the offset velocity range [−2.500; +1.500] km s −1 (blue triangle). We compare their Lyman-α luminosity-weighted number densities with the expected values for field galaxies at comparable redshift according to the Lyman-α luminosity functions measured with MUSE in wide 27 (gray line and yellow shaded region, representing the ±1σ uncertainty on the fitting parameters) and deep 26 surveys (red points and black lines, representing the ±1σ errors). The number density of Lyman-α-emitting galaxies, n LyE , associated with W0410-0913 is a factor of about 14 higher than the value (<n LyE >) expected for field galaxies with comparable Lyman-α luminosity. The green and blue error bars represent the range of Lyman-α luminosities covered by companion galaxies in the whole and denser volume, respectively.
ALMA observations of the internal ISM structure. W0410-0913 has been extensively observed by ALMA targeting several rest-frame FIR lines in different bands. While existing ALMA observations cannot provide information on the companion galaxies, they succeed in describing several properties on scales of the ISM, complementing the information retrieved by MUSE on the CGM environment. Studying the emission from the rotational transition of carbon monoxide CO , observed in Band-4 with a high resolution of about 0.2″ (corresponding to 1.4 kpc at this redshift), we find a smooth velocity gradient (see a map of the line-of-sight velocity in Fig. 4a) with velocities ranging  between approximately ±500 km s −1 . We perform a kinematic model of W0410-0913 using the 3D tilted ring model fitting code 3D Barolo (see details on the procedure in Methods) and obtain models of the rotational velocity (v rot ) and velocity dispersion (σ) radial profiles across the whole CO(6-5) emitting region (Fig. 4a). The v rot and σ radial profiles are shown in Fig. 4b. We find that, after a declining trend in the inner 1−2 kpc, v rot (r) follows a flat curve (within the uncertainties) reaching a value of about 500 km s −1 at a radius of~6 kpc. The σ(r) profile instead peaks at the center (about 200 km s −1 in the inner 1−2 kpc) and slowly declines outward, exhibiting high values of σ ≳ 100 km s −1 up to about 5−6 kpc. In the lower panel of Fig. 4b we show, as a function of the radius, the ratio between rotational velocity and velocity dispersion, v rot /σ, which is typically used as an indicator of whether a system can be gravitationally supported by rotation or by turbulence 39 . We find that, at each radius, the rotational velocity is higher than the dispersion velocity, with an average (median) value of v rot /σ = 3.5 ± 1.8 (v rot /σ = 3.0 ± 0.5). While W0410-0913 appears to be more turbulent than a few recently reported extreme dynamically cold disks at z > 3 [39][40][41] (that show a surprisingly large v rot /σ~10), the v rot /σ values that we estimate are consistent with state-of-art galaxy formation models 42,43 , indicating that the molecular gas emission observed in our target arises from a massive rotationally supported fast rotating disk.

Discussion
ALMA observations of W0410-0913, zooming into the kpc-scale properties of its ISM, interpreted in the context of the CGM-scales studied with MUSE, challenge the scenario, suggested by models, that repeated galaxy encounters in the dense environment might have been efficient in building up the central galaxy. Indeed, multiple mergers in luminous systems are often shown to have a disruptive impact on the kinematics of the ISM 17,35,44 , which, in the case of W0410-0913, appears to be settled into a smooth rotation-dominated structure.
We speculate that a possible mechanism capable of solving this tension is that the mass growth of W0410-0913 could be sustained by minor mergers with relatively less massive companions, that do not significantly perturb the disk-like structure 43,45 . In this case, the accretion of discrete small satellites might have played a role in originating the observed high-velocity dispersion in W0410-0913 through gravitational instabilities. Also, minor mergers are usually invoked by models as an efficient channel of bulge formation at the core of massive disk galaxies 46 . In the case of W0410-0913, the inner peak (<1−2 kpc) of the CO(6-5) v rot (r) and the high central concentration of flux in the intensity map could be interpreted as a tracer of a bulge-like structure, that minor mergers might have contributed to form. However, we note the existence of an alternative (but not exclusive) scenario in which the energy injected by the prominent AGN activity into the surrounding gas might contribute to the observed turbulence throughout the disk, as suggested by similar observations of turbulent ISM revealed through spectroscopic imaging of molecular gas in another distant hyper-luminous obscured quasar 13,47 .
Altogether, the multi-wavelength and multi-scale observations presented here, obtained with some of the most powerful telescopes available (including the first observations targeting a Hot DOG with MUSE) provide key information on the cosmic assembly of a distant hyper-luminous hot-dust obscured galaxy. Data support that the mass build-up of W0410-0913 takes place in a dense circumgalactic environment, but we cannot identify robust signs of merger that could sustain the growth of the massive galaxy. Instead, we find that gravitational encounters, if any, can have only a modest impact on shaping the morpho-kinematic properties of the ISM. Our results serve as a benchmark to calibrate models of the rare and short-lived populations of Hot DOGs and progress in the understanding of their evolutionary transition towards luminous quasars and giant galaxies at cosmic noon.

MUSE observations
We have processed the individual exposures from both datasets with standard data reduction techniques using the ESO-MUSE Pipeline Version 2.8.4 49 . In detail, for each exposure, we used the MUSE_BIAS, MUSE_FLAT, and MUSE_WAVECAL pipeline recipes 49 to perform bias subtraction, flat fielding, twilight and illumination correction, and wavelength calibration. In order to ensure an accurate relative astronomy, we registered the individual exposures using the position of point sources in the field. This step is needed since the spatial shifts introduced by the derotator wobble 50 can cause offsets of a few arcsecs.
We performed the next stages of data reduction using the tools provided by the CubExtractor software package 37 (Version 1.8), which has been specifically developed to improve the data quality for the detection of diffuse, low-SB emission in MUSE datacubes 22,23,51,52 . In detail, we improved the flat fielding by using the tool CubeFix 37 , which performs an additional correction to homogenize the illumination across the field and as a function of wavelength. We performed the sky subtraction using CubeSharp 37 , which can reach a better accuracy than the MUSE pipeline 22 . We then combined the corrected and skysubtracted individual exposures from both datasets using an average 3σ-clipping. The final datacube samples the instrument FOV of 1 0 × 1 0 in pixels of size 0.2″. Each pixel contains a spectrum covering the wavelength domain 4750−9350 Å, and individual channels have a spectral resolution of 1.25 Å.
Our combined MUSE observations, comprising a total telescope time of 2h, reach a (2σ) SB limit of 5 × 10 À19 erg s À1 cm À2 arcsec À2 in 1 arcsec 2 aperture and in a single spectral channel (1.25 Å), at wavelengths adjacent to 5632 Å, corresponding to the expected Lyman-α in W0410-0913 (systemic redshift of z = 3.631; see next Sections). For reference, the SB limit on a pseudo-narrow band (NB) image of 30 Å (a common size for filters used in past NB surveys), obtained by collapsing 24 wavelength layers close to the location of Lyman-α, is 2:8 × 10 À18 erg s À1 cm À2 arcsec À2 in 1 arcsec 2 (2σ).
At the systemic redshift of W0410-0913, z = 3.631, traced by its molecular gas emission, we do not find any Lyman-α emission. The average flux, computed within a spectral window of Δv = 1500 km s −1 centered around the expected line center at 5632 Å, is consistent with zero within the uncertainties,~3(±30) × 10 −20 erg s −1 cm −2 Å −1 (see Fig. 5b), and by producing pseudo-NB images at the expected position of Lyman-α (collapsing the continuum-subtracted datacube along the wavelength axis using different spectral widths, from 5 to 30 Å) we do not find any sign of emission. However, a narrow emission line (FWHM~400(±80) km s −1 ), emitted at the same location of W0410-0913, is detected in the spectrum with a blueshift of about −1500(±150) km s −1 with respect to z = 3.631 (see Fig. 5). We interpret this line as a Lyman-α emission from W0410-0913, emitted at different velocities than the systemic redshift traced by the molecular gas. This interpretation is supported by the blueshift of about −1400 km s −1 (with respect to the systemic redshift adopted in this work) of the narrow component of the non-resonant line [OIII] λ 5007 from a rest-frame optical spectrum of W0410-0913 observed with KECK/NIRES 5 . These findings suggest that the rest-frame UV and optical narrow emission lines in W0410-0913 are emitted at different velocities than the molecular gas. An alternative, but less likely (because of the reasons discussed above) interpretation is that the narrow rest-frame UV emission line that we find in the MUSE spectrum traces Lyman-α emission produced by a nearby companion galaxy that is part of the overdense environment (see subsection The circumgalactic environment of W0410-0913 observed with MUSE), and is located along the line of sight of W0410-0913.
We compute the spatial extension of such emission using the tool CubEx (from the CubExtractor package 37 ). CubEx is extensively used to perform automatic detection and extraction of extended sources in MUSE datacubes 22,23 , and it is based on a 3D extension of the connected-labeling-component algorithm with union finding, usually used in classical binary image analysis 22 . We feed CubEx with a subcube that we extract to have a wavelength range of Δλ~100 Å centered around the peak of Lyman-α emission (4388 Å), and from which we remove continuum sources by adopting a median-filtering approach. We also apply a spatial Gaussian filtering of 2 pixels (0.4″) and a spectral smoothing with a Gaussian filter size of 1 layer (1.25 Å) to bring out extended but narrow features. Running CubEx with a SNR threshold for detection of SNR > 3 in individual datacube elements (voxels) we detect a nebula consisting of approximately 1600 connected voxels and a maximum spectral width, as defined in the 3Dsegmentation-mask, of Δλ~15 Å. From the optimally extracted image of the nebula, derived by applying the 3D-segmentation-mask to the datacube (see subsection A study of the 3D distribution of gas around W0410-0913, for details on the advantage of this method), we measure a spatial extension, defined as the maximum projected linear size, of about 30 kpc (see Fig. 6). This indicates that, under the reasonable assumption that such emission could be considered a Lyman-α nebula associated with the Hot DOG, it would be smaller than several nebulae

The circumgalactic environment of W0410-0913 observed with MUSE
To study the environment of W0410-0913 over an area of about 450 kpc × 450 kpc (corresponding to the MUSE 1 0 × 1 0 FOV), we search for Lyman-α-emitting companion galaxies in a continuum-subtracted sub-cube (produced through median-filtering) of about 380 Å centered around the expected Lyman-α emission from W0410-0913, corresponding to offset velocities of [−10.000; + 10.000] km s −1 with respect to its systemic redshift.
Our method can be summarized as follows. We inspect the whole wavelength domain by producing mini-cubes with spectral widths ranging from a minimum allowed number of 2 channels (2.5 Å; to avoid single-channel spikes of noise) to a maximum of 20 channels (about 25 Å; to allow for~3 × typical FWHMs of Lyα-emitters at this redshift). We collapse each mini-cube along the wavelength direction to obtain pseudo-NB images that we scan through their whole spatial domain by measuring the aperture photometry of 1″−sized regions centered around each pixel. We consider as candidate detections the regions with a significance of their measured flux of SNR > 3 (where the noise is extracted from the local background), we extract a spectrum, and we save their position (RA and DEC of the central pixel), central wavelength and spectral width. For multiple candidate detections arising from the same physical source, namely two or more outputs found at about the same location (within 0.5″, corresponding to about half of the size of the PSF) and wavelength (within 4 Å), we keep the candidate that has higher SNR in the spectrum. We then perform a visual inspection of the pseudo-NB images and spectra of the candidate detections to make sure that there is no contamination by sky residuals. Also, for each candidate detection, we visually inspect the original non-continuum subtracted full datacube (covering the whole MUSE spectral domain) to check whether the candidate detections show the clear presence of multiple (at least two) emission lines at wavelengths not included in the sub-cube used for studying the environment of W0410-0913. In these cases, the sources would be line emitters at lower redshifts (emitting e.g., [ and Hα detected in the full spectrum). We assume that the remaining sources, showing a single emission line in the sub-cube centered around the redshift of W0410-0913 are Lyman-α-emitting galaxies in the surroundings of W0410-0913.
We find 24 Lyman-α emitters within a box with a diameter of about 400 kpc centered around W0410-0913. This result is confirmed by an independent analysis that we perform using CubEx. Although CubEx is usually exploited to extract diffuse extended emission, it can be used to mimic a line finder by decreasing the threshold of the minimum allowed number of connected voxels and by increasing the threshold of SNR for the detection of individual voxels. In detail, by adopting a SNR > 4 threshold over a lightly smoothed cube (0.2″ pixels in the spatial domain, and no filtering on the wavelength direction), we find the same number of sources of the method described above, with only tiny differences (≲0.3″) in the spatial location of the peak of a few objects, likely due to the different pre-processing of the cube.
In Supplementary Figs. 1 and 2, we show the spectra of these sources (extracted from a region with a radius of 1″) centered around the central wavelength of their Lyman-α, that we compute through a non-linear least-squares Gaussian fit. A selection of these spectra is shown in Fig. 7a. From the modeled Lyman-α peaks we measure the offset velocities with respect to the redshift of W0410-0913, and find values in the range [−6.000; + 10.000] km s −1 , spatially distributed in the MUSE FOV as shown in Fig. 1b. We note that radiative transfer effects are known to shift the peak of Lyman-α emission in the spectra of star-forming galaxies even by a few hundreds of km s −1 compared to the systemic velocity 56 , though this effect is unlikely to have a dominant impact on the suggested interpretation of the observations, since we measure offset velocities up to thousands of km s −1 .
We also produce pseudo-NB images of each object by collapsing the datacube across the wavelength direction within spectral ranges defined by [λ cen − FWHM:λ cen + FWHM], where the central wavelength λ cen and FWHM are defined by the spectra shown in Supplementary  Figs. 1 and 2. A composition of the derived pseudo-NB images, showing the projected spatial distribution of the Lyman-α-emitting galaxies around W0410-0913, is shown in Fig. 1a.
As shown in Fig. 1c, the distribution of offset velocities is centered around W0410-0913, confirming that the massive luminous galaxy is located at the dynamic center of its dense environment and dominates its gravitational potential. Most of the companions (18 out of 24) show offset velocities in a smaller range of about [−2.500; +1.500] km s −1 .
With the caveat of partial information provided by the measurement of projected line-of-sight velocities, it is likely that these galaxies are more closely interacting with the central massive object.
We dub the Lyman-α-emitters using alphabetically ordered letters from #A to #X, according to their velocity offsets, from the largest negative to the largest positive. We compute the Lyman-α fluxes of the companion galaxies by integrating the Gaussian model of their spectra, and convert them to Lyman-α luminosities (L Lyα ), obtaining values in the range logðL Lyα =½erg s À1 Þ~41. 8-42.65. A summary of the measured properties of the Lyman-α-emitters is reported in Table 1.

A statistical assessment of the overdensity
To provide a statistical assessment of our finding, we perform the same analysis discussed in the previous subsection on other four 380 Åsized sub-cubes adjacent to the sub-cube used to study the environment of W0410-0913, in the large wavelength range 4750−6500 Å (corresponding to offset velocities of [−50.000; +50.000] km s −1 ). This wavelength range is optimal because the sensitivity of the instrument varies only by a factor less than 25%, and it is not strongly affected by sky lines emission (see a description of MUSE performances at https:// www.eso.org/sci/facilities/paranal/instruments/muse/inst.html). After checking for contamination from sky-residuals and multiple line emitters at intermediate redshifts (the number of the latter is found to be between 0 and 2 per sub-cube), we find that the number of detected candidate Lyman-α emitters is either 1 and 2, which confirms the overdense nature of the environment in the surrounding of the Hot DOG.
As an additional sanity check, using the same method, we scan through the negative version (obtaining by changing the sign to each pixel) of the sub-cubes defined above. In this case, the visual inspection following the line finder algorithm is limited to a check of the contamination by sky-residuals. We find that the number of detections is 0 in all the sub-cubes. This excludes the possibility that a fraction of the positive sources could be resulting from spikes of noise and thus be spurious.

The physical significance of the overdensity
To quantify the significance of the discovered overdensity we compute the Lyman-α luminosity-weighted number density (Φ) of Lyman-αemitting galaxies surrounding W0410-0913 and compare it with the value predicted by the most updated Lyman-α luminosity functions measured for field galaxies in wide and deep surveys with MUSE 26,27 .
The 1 0 × 1 0 projected area (i.e., the MUSE FOV) over which we find the companion galaxies corresponds to a comoving area of 4.2 Mpc 2 . In the whole volume case, the redshift-side has a comoving size of about 190Mpc, which yields to a cosmological comoving volume of 800 Mpc 3 and a number density of Lyman-α emitters of n LyE3 × 10 −2 Mpc −3 . In the dense volume case, the redshift-side has a smaller comoving size of 50Mpc, yielding to a cosmological comoving volume of about 210 Mpc 3 and a number density of Lyman-α emitters of n LyE8 .5 × 10 −2 Mpc −3 . To convert n LyE into Φ we divide by the range of  Table 1). Lyman-α spectra, extracted from circular regions with a radius of 1″ and centered around the Lyman-α peaks, are shown in gray. The wavelength range in the x-axis has a width of 120 Å and is centered around the central wavelength obtained through a Gaussian fit (Gaussian models are shown in blue). Legends report the redshift of each source (derived through Lyman-α), their FWHM, and their offset velocity with respect to the systemic redshift of the central luminous object. Gray shaded areas represent the regions affected by sky lines. A compilation of MUSE spectra for all the Lyman-α-emitting galaxies around W0410-0913 is shown in Supplementary Figs. 1 and 2. b HST cutouts (4″ × 4″) of the Lyman-α shown above, detected with the WFC3/F160W filter (colored images). The 2σ contours of pseudo-NB Lyman-α images obtained with MUSE are overlaid in black. Red dashed circles represent the circular apertures with 1″-sized radii that we use to extract the HST flux and estimate the SFR. The cutouts for all the HST-detected companion galaxies are shown in Supplementary Fig. 3.
Lyman-α luminosities of galaxies in both volumes, namely Δ logðL Lyα =L Þ = 0:86 for the whole volume and Δ logðL Lyα =L Þ = 0:75 for the dense volume. We obtain Lyman-α luminosity-weighted number densities of Φ whole = 3:5 × 10 À2 dlogðL Lyα =L Þ À1 Mpc À3 and Φ dense = 1:1 × 10 À1 dlogðL Lyα =L Þ À1 Mpc À3 . These values are higher than the expected values for field galaxies (see Fig. 2). In the following we focus on the dense volume, which is more likely to contain companion galaxies gravitationally associated with W0410-0913. In detail, at the median Lyman-α luminosity of galaxies in the dense volume, logðL Lyα =L Þ = 42:2, the MUSE-wide survey 27 , which observed Lyman-α-emitters over a large area of >20 arcmin 2 , yields an expected Φ field = 10 Àð2:1 + 0:37 À0:47 Þ dlogðL Lyα =L Þ À1 Mpc À3 , where the uncertainty represents the errors on the best-fitting Schechter function parameters 27 . Comparing this value with Φ dense , we find that the circumgalactic environment around W0410-0913 has a denser concentration of Lyman-α-emitting galaxies than the average volume, with an estimated overdensity factor of ρ CGM = Φ dense Φ field = 14 + 16 À8 . We note that this calculation is conservative since it is based on the assumption that offset velocities with respect to the central object are tracing relative distances, while it is more likely that the observed dispersion of offset velocities (at least for objects in the dense volume) is affected by gravitational attraction (see Article and previous subsections). Thus the effective values of Φ in our field are likely to be higher, e.g., up to a factor of about 20 times higher in the dense volume case under the extreme assumption that the 18 Lymanα emitters are contained in a cube with a side of 60″ centered around W0410-0913.
Despite its high value, we find that the measured overdensity of Lyman-α emitters in the dense volume, ρ CGM = 14 + 16 À8 , is consistent with cosmological predictions in the proximity of extreme, luminous objects like W0410-0913.
To check this, we first estimate how rare the fluctuation of matter density associated with ρ CGM would be. The first step is to consider that in general the spatial clustering of observable galaxies does not precisely mirror the clustering of the bulk of matter in the Universe. This caveat is often addressed through the definition of a galaxy clustering bias, which describes the relationship between the spatial distribution of galaxies and the underlying dark matter density field 57 . Thus, by using a standard assumption of linear bias 58 , we convert the observed overdensity of Lyman-α emitters ρ CGM into an estimate of the underlying matter overdensity, ρ M , by dividing it by the clustering bias parameter of Lyman-α emitters, b. By using an updated value of b = 2.8 ± 0.4, recently computed through a 3D clustering analysis of hundreds of Lyman-α emitters at z ≲ 4 in the large MUSE-Wide survey 59 , and consistent with previous similar studies 58,60 , we obtain an estimate of the matter overdensity in the surroundings of W0410-0913 of ρ M = 5 + 7 À3 . To estimate the rarity of such a fluctuation, we start from the computation of the linear rms density fluctuation 61 , σ LN , predicted by the standard ΛCDM cosmology. For the dense volume over which we compute ρ M , i.e., about 210 Mpc 3 , we find a value of σ LN~1 :6. The next step is to compute the probability associated with observing an overdensity ρ M given a certain density fluctuation. On the physical scales of our system, the density field is non-linear, and a possible approximation for the true density distribution P(ρ) is a log-normal model 62 , for which the σ LN~1 :6 can be interpreted as the rms in ðρÞ, with a mean in this quantity of À  In the first, second, and third columns we report the IDs and the coordinates (Right Ascension and Declination) associated to the companion galaxies detected through Lyman-α in the surroundings of W0410-0913. The coordinates correspond to the centroid of Lyman-α emission in the pseudo-NB images. Redshifts, FWHMs, and LLyα (in the fourth, fifth and sixth column, respectively) are derived from the MUSE spectra as described in the Methods, and their errors represent the uncertainties provided by the Gaussian models. The SNR values in the seventh column are computed as the ratio between the integrated flux in the MUSE spectra and the noise associated to the width of the line (for each source, we considered spectral regions as large as 2 × FWHM). In the last column, the SFRs and their uncertainties are based on the HST images and are computed as discussed in the subsection Measure of the SFR from HST in Methods. P(ρ M ), we find that the comoving density of objects like the central galaxy would need to be about n halo~P ðρ M Þ 210 Mpc 3~1 0 À3:7 ± 0:6 Mpc À3 . This values is consistent, within the uncertainties, to the predicted comoving density of halos hosting central massive galaxies like W0410-0913, with maximum velocities of v max = 500 ± 70 km s À1 (see the radial profile of the rotational velocity in Fig. 4b) at z = 3.631, which we find to be n halo~1 0 −3.8±0.5 Mpc −3 , as we compute integrating the halo mass function provided by the large N-body simulation MDPL2 63 , included in the open-access CosmoSim database (https://www. cosmosim.org/cms/simulations/mdpl2/).

Measure of the SFR from HST
We estimate the SFR of the Lyman-α-emitting companion galaxies around W0410-0913 using archival HST data, that we retrieve from the public Barbara A. Mikulski Archive for Space Telescopes (MAST) portal (https://archive.stsci.edu/). A 2 0 × 2 0 field around our object was observed for a total exposure time of 2400 s in October 2012 during the HST Cycle 20 (Program ID: 12930) exploiting the filter F160W on the Wide Field Camera 3 (WFC3), centered at about 15436 Å, which corresponds to~0.3 μm at the redshift of interest. We extract the fluxes in the HST WFC3/F160W image from circular apertures with 1″-sized diameters centered around the peak emission, that we search within a radius of 0.5″ from the Lyman-α peak (obtained through a 2D Gaussian fit of the pseudo-NB images; see Fig. 1a), to allow for slight spatial offsets between Lyman-α and the rest-frame UV/ optical continuum 64 . We then derive the luminosity density L ν and convert it in SFR by using a standard rest-frame UV/optical L ν −SFR calibration and assuming a Chabrier initial mass function 65 . Although the available HST data are shallow (with a 3σ flux limit corresponding to a SFR 3σ~9 M ⊙ yr −1 ), we detect 10 galaxies at SNR > 4, with star formation rates ranging between SFR~12−100M ⊙ yr −1 .
In Supplementary Fig. 3, we show the HST cutouts for the companion galaxies detected with the WFC3/F160W filter, along with the Lyman-α contours from MUSE. A selection of HST maps is shown in Fig. 7b. The SFRs are reported in Table 1.
To obtain a better constraint of the average star formation rate in the 14 companion galaxies not detected by HST we perform a median stacking of their 5″ × 5″ cutouts centered at the Lyman-α peaks. Prior to stacking, we run a sigma-clipping algorithm (with a 3σ threshold on individual pixels) over each cutout to mask possible continuum sources outside the central circular area of 1″ diameter. This is necessary to guarantee an uncontaminated estimate of the noise level in the final image, that would be otherwise affected by bright field objects. The noise level in the stacked image decreases down to a corresponding 3σ SFR limit of SFR 3σ~2 .5M ⊙ yr −1 , enabling us to measure with a SNR~3.6 an average SFR~3M ⊙ yr −1 , suggesting that a non-negligible star formation is taking place in the companion galaxies not detected by HST too.
A study of the 3D distribution of gas around W0410-0913 Extended diffuse Lyman-α emission, mainly produced through recombination radiation following photoionization (often referred to as fluorescence) powered by UV-bright sources 24 , is a commonly used tracer to observe the cool gas (with temperatures of approximately T = 10 4 K) in the CGM around galaxies. At z ≳ 2−3, where the Lyman-α is redshifted into the optical domain, great advances have been recently made with wide-field integral field spectrographs at large telescopes (like MUSE, or the Keck Wide-Field Imager), which have identified ubiquitous Lyman-α nebulae around quasars, extending up to hundreds of kpc 22,23 , and Lyman-α halos around star-forming galaxies, extending up to tens of kpc, between factors of 5−20 larger than the rest-frame UV sizes 66 .
We search for extended Lyman-α emission around the companion galaxies of W0410-0913 to study the morphological and kinematic distribution of gas in the circumgalactic environment of W0410-0913. We run CubEx over the sub-cube in which the 24 Lyman-α-emitters are found, using the same pre-processing and setup of parameters adopted for analyzing the Lyman-α nebula (see previous subsections), namely a SNR > 3 threshold for voxel-detection in a smoothed cube with a Gaussian spatial filter of 0.4″ and a Gaussian spectral filter of 1.25 Å. This setup is different from that exploited in the line-emitters searches, and it is tailored at extracting faint and narrow features. We produce optimally extracted Lyman-α flux maps of the sources by selecting the voxels detected in their 3D-segmentation-masks and integrating their fluxes along the wavelength axis. This technique is operationally equivalent to combining different pseudo-NB images in which the size of the spectral filter is optimized for each pixel. The advantage of this approach, enabled by integral field spectrographs, is that one can maximize the SNR of the whole emitting region by revealing both kinematically narrow and broad features of a source at the same time. In contrast, imaging through classical NB filters with a fixed width would produce either a noise-dilution or an underestimation of the signal emitted.
A collection of the optimally extracted images of the companion galaxies is shown in Fig. 3a. We find that some of the luminous companion galaxies in the denser dynamic volume around W0410-0913, with offset velocities in the range [−2.500; +1.500] km s −1 , show extended diffuse Lyman-α envelopes (not revealed, or barely detected, in the pseudo-NB images; see Fig. 1a) that are suggestive of connecting filamentary structures on scales up to ≳100 kpc (see Fig. 3). We use the 3D-segmentation-masks provided by CubEx to produce also the 2D maps of the first moment (moment-1) of the flux distribution in its spectral domain. These maps give an indication of the centroid offset velocities with respect to the systemic redshift of W0410-0913 at each spatial location. A collage of such kinematics maps is showed in Fig. 3b.

ALMA data
In this work, we use archival ALMA observations of W0410-0913 targeting the emission-line spectra arising from rotational transitions of carbon monoxide CO(4-3), CO , and CO , and the atomic car- The pipeline calibration outputs did not show any issue, and thus we did not perform any additional flagging. We obtained continuumsubtracted visibilities using the CASA task uvcontsub, and we generated the continuum-subtracted datacubes of the spectral windows containing the emission lines by running the CASA cleaning algorithm clean with 500 iterations and a SNR threshold of 3. In order to maximize the sensitivity, we adopted a natural weighting of the visibilities for all the three datasets. The resulting synthesized beam in the CO(4-3) datacube is 0.7″ × 0.5″ (position angle, PA = 83°), 0.22″ × 0.18″ (PA = −69°) for CO , and 0.6″ × 0.4″ (PA = 80°) in the datacube containing CO(7-6) and [CI]. We adopted pixel sizes of 0.1″, 0.03″, and 0.1″, and spectral bins of 50 km s −1 , 60 km s −1 , and 50 km s −1 (which enable us to have both a good sampling of the lines and enough SNR per individual channel) for CO(4-3), CO(6-5) and CO(7-6) + [CI] respectively. The average sensitivities reached in the cubes (in spectral regions close to the line frequencies) are 0.7 mJy beam −1 in a spectral channel of 50 km s −1 for CO(4-3), 0.08 mJy beam −1 in a spectral channel of 60 km s −1 for CO  and 0.8 mJy beam −1 in a spectral channel of 50 km s −1 for CO(7-6) + [CI].
In Fig. 8 we show the CO(4-3), CO , CO(7-6) and [CI] spectra of W0410-0913. These spectra are extracted from a spatial region defined by the connected pixels with SNR > 2 in the moment-0 maps, that we first obtain by collapsing the cube spectral channels in the range [f cen −FWHM: f cen + FWHM], where f cen (the central frequency) and FWHM are derived by a Gaussian fit of a spectrum extracted from a 1″sized circular aperture centered at the expected position of the object. Such a 2-iterations procedure helps in case of resolved observations to ensure that the whole line-emitting region is accounted for in the spectrum while avoiding signal dilution. Through a non-linear leastsquares Gaussian fit of the CO(6-5) line, which has the highest SNR among the available data, we derive a systemic redshift of z CO(6−5) = 3.631 ± 0.001. This value is consistent with an independent measurement of z CO(6−5) recently reported 18 (z = 3.6301), and the tiny difference is probably due to a slightly different definition of the extracting aperture. The velocity axes of the spectra in Fig. 8 are aligned to the common systemic redshift of z CO(6−5) = 3.631. From our Gaussian fit we derive large FWHMs of about 800 km s −1 (with average associated uncertainties of 70 km s −1 ) in CO(4-3), CO  and CO .
[CI] appears to be narrower (FWHM~600 ± 90km s −1 ), although this might be an effect of spectral blending with the more luminous CO . Broad molecular line profiles in W0410-0913 were also reported in previously published ALMA observations of CO(4-3) 19 , and [CII] 18 , although in both cases the instrument spectral setup was tuned on the optical redshift 2 , which is blueshifted with respect to the systemic redshift as traced by the molecular gas, and the lines, close to the edge of the spectral window, were only partially detected. Such broad emission lines in the integrated spectra extracted from the whole galaxy reflect the rapidly spinning nature of its disk (see subsection Kinematic modeling of the CO(6-5) emission in the ISM of W0410-0913).

Study of the CGM environment with ALMA
Unfortunately, none of these datasets, designed to observe the cold gas emission from the hyper-luminous W0410-0913, provides information about the rest-frame FIR line properties of the companion Lyman-α emitters that we discover with MUSE in its CGM. However this is not surprising, and it is mainly due to a combination of insufficient sensitivity and too high resolution (that dilutes the flux over several beams), as we test by running simulations of the estimated telescope time on the ALMA Observing Tool (OT; version released for Cycle 8). In detail, we first estimate the total molecular gas mass (M mol ) from the SFR measured with HST by using standard M mol -SFR scaling relations calibrated at high-z 68 . We then convert M mol into total CO fluxes by adopting a conservative Milky-Way-like conversion factor between M mol and the CO(1-0) luminosity, α MW CO = 4:3, and standard Milky-Waylike intensity ratios between CO(1-0) and the rotational transitions observed by ALMA 69 . We note that the adopted value of α CO has been calibrated at solar metallicities, and should be considered conservative in our case since α CO typically increases at lower metallicities 70 , and high-z galaxies tend to be less metal-rich at a fixed stellar mass 71 .
Assuming projected sizes of about 1″ and line-widths of 300 km s −1 , we find that the estimated total ALMA time to observe the integrated CO spectrum in a galaxy with SFR = 30 M ⊙ yr −1 (the median SFR of the HSTdetected companions), with at least 3σ per channel in at least 3 spectral bins, at the same angular resolution of the available observations, would be about 1.6h for CO(4-3) and CO , and about 24.5h for CO , between a factor of 6−8 longer than the observed datasets. Other additional limitations derive from the ALMA smaller fields of view (the primary beams range from about 30″ in Band 5 to about 45″ in Band 3), and the narrower spectral coverage (sidebands ranging from 5000 km s −1 in Band 5 to 9000 km s −1 in Band 3) with respect to MUSE data. The combination of these effects reduces the number of companion galaxies that could be possibly detected (up to 60%, in the best case) and, together with the sensitivity issue discussed above, limits the information that can be derived through stacking analyses. Indeed we find signals that are consistent with noise in the stacked images that we produce, for each dataset, by combining moment-0 maps centered at the expected positions of the galaxies (using the peak of Lyman-α emission as a prior), obtained by collapsing three different velocity ranges, 300, 400, 500 km s −1 , centered around the expected systemic frequencies (using the redshifts derived from the Lyman-α as a prior).
Kinematic modeling of the CO(6-5) emission in the ISM of W0410-0913 The deep and highly resolved CO(6-5) observations of W0410-0913 are best suited to study the molecular gas kinematics in W0410-0913, that Article https://doi.org/10.1038/s41467-022-32297-x we model using the 3D tilted ring model fitting code 3D Barolo 72 (3D-Based Analysis of Rotating Objects from Line Observations). This technique has the advantage that it directly models the 3D cube taking beam-smearing effects into account and it has been extensively applied to observations of high-redshift (z > 2) galaxies 73 . The fit routine is performed on the continuum-subtracted CO(6-5) datacube and can be summarized as follows.
First, we isolate the signal by producing a 3D signal mask of the emission line through the 3D Barolo SEARCH tool (based on DUCHAMP 74 ). In detail, this algorithm identifies pixels in the datacube with intensities above a certain SNR (SNR upper = 3.5, in our case), and then it searches around these signal peaks for emission above a second SNR threshold (i.e., SNR lower = 2.0, in our case). The 3D Barolo best-fit results are highly dependent on the overall geometry of the model (i.e., number and width of the rings) and on the initial estimates of the parameters. Thus, once the signal is isolated we create the moment maps needed for the initial guesses. Using the CASA toolkit task image.moments we collapse the cube over the spectral channels containing the line and produce a moment-0 map, which we then fit with a 2D Gaussian (using the CASA task image.fitcomponents) to obtain the FWHM of the major and minor axes, the position angle and the central position. We also produce a velocity dispersion map (i.e., moment-2) including only the pixels that belong to the 3D mask identified before. Then, we fit a tilted ring model to the line emission by using the 3D Barolo key function 3DFIT. We adopt a value of twice the FWHM of the major axis of the restoring beam as a maximum model  In panel a, b, c, we show the moment-0, moment-1, and moment 2 of the CO(6-5) data. For each of these panels, the three columns denote (from left to right) the observed data, models (created with 3D Barolo), and the corresponding residuals. The black solid line in the lower-left corner of the moment-0 panel shows a 5 kpc physical scale. The solid lines in the second row represent the kinematic major axis, while the dashed lines represent the minor axis. The black contours in the moment-0 maps are defined as in Fig. 4. In panel d we show the major axis (left) and the minor axis (right) PVDs. For each PVD, the data is shown by the background color, while the overlaid contours represent the model at 20%, 50%, and 80% of its maximum value. The ALMA beam, colored in red, is shown in the lower-left corners. radius. As the value for the width of each ring, we adopt the FWHM of the minor axis, in order to not over-resolve the data 73 . This corresponds to a number of 12 model rings. The model central position is fixed to be coincident with the centroid of a 2D Gaussian fit of the moment-0 map, and we use the maximum value in the moment-2 map as the initial guess of the velocity dispersion. A thin disk, with a scale height of 0.01″, is assumed. We adopt an initial guess for the inclination of 45°, and we allow it to vary between 10 and 80°. Starting from this set of guesses, 3D Barolo produces a model of the first ring by randomly populating the volume with discrete emitting gas clouds in order to reproduce the position angle, the rotational velocity, etc of the initial estimates. The resulting model datacube is then convolved with the ALMA synthesized beam and rescaled so that the velocity-integrated intensity of each spaxel in the data and model cubes are equal. Then, the values of the pixels in the data and model datacubes are compared, and each parameter is varied until the absolute residual (i.e., |model-observation|) is minimized. When this process reaches a minimum, the next ring is analyzed. During the first run 3DFIT provides best-fit results for the rotational velocity, velocity dispersion, systemic velocity, inclination, and position angle. The three latter variables are averaged across all the rings, and the whole fit routine is run again, with only the rotational velocity and velocity dispersion as variables. Among the final 3D Barolo outputs, we obtain morphological parameters (e.g., central position), and kinematic parameters such as inclination, position angle, velocity dispersion, and rotation.
In Fig. 9 we show that our model is able to reproduce the kinematic properties of the target to an excellent level, with small residuals in the moment maps and a good agreement between the model and data position-velocity diagram (PVDs). The galaxy has a well-defined rotation, with rotational velocities that range between v rot~± 500 km s −1 , as can be seen from the moment-1 map, showing the spatial distribution of line-of-sight velocities, and from the PVDs extracted along the major and minor axes (see Fig. 9). Despite the fast coherent rotation, the velocity dispersion is still large, in analogy with other recent ALMA observations of Hot DOGs reporting turbulent ISM 18 . Nevertheless, as discussed in the Article, we find that v rot /σ > 1 at each ring in our model (see Fig. 4). We find morphological and kinematic position angles of 133 ± 5°and 133.8 ± 4.9°respectively, and morphological and kinematic inclination on the plane of the sky of 49.5 ± 3.3°and 50.1 ± 13.1°.

Code availability
The MUSE pipeline recipes, used to reduce the MUSE data, are publicly available and the software distribution can be obtained at https://www. eso.org/sci/software/pipelines/muse/. The software CASA, used to process the ALMA data, is publicly available and can be downloaded at https://almascience.eso.org/processing/science-pipeline. The code used to model the molecular gas kinematics, 3D Barolo, is publicly accessible through github (https://editeodoro.github.io/Bbarolo/), with a full documentation hosted at https://bbarolo.readthedocs.io/en/ latest/. The package used for the 3D extraction of faint diffuse emission in the MUSE datacube, CubExtractor, is not publicly available, though the reader interested in using this code can contact Sebastiano Cantalupo at the mail address sebastiano.cantalupo@unimib.it.